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I .  INTRODUCTION 


The  application  of  microprocessors,  which  began  early 
in  1972,  is  rapidly  becoming  a  major  computational  activity. 
The  control  of  traffic  lights ,  smart  terminals  and  cash 
registers,  printer  controls,  etc.,  were  some  of  the  first 
applications.  Future  applications  are  almost  unlimited, 
starting  with  present  implementations  of  complex  plant 
process  controls  and  numerically  controlled  machines. 

Because  of  their  small  size,  light  weight,  low  power 
requirement  and,  most  important  of  all,  low  cost,  micro¬ 
computer  systems  are  ideally  suited  for  use  on  board 
submarines.  Data  processing  with  the  aid  of  microcomputers 
is  possible  in  almost  all  technical  and  organizational 
fields  in  a  submarine,  such  as: 

1.  sensors  and  signal  processors, 

2.  fire  control  and  weapons  control, 

3.  communications  and  navigation, 

4.  central  command,  information  display  and  administration, 

5.  ship  handling  and  automatic  control. 

This  thesis  developed  an  algorithm  for  controlling  depth 
and  pitch  of  a  submarine  using  a  microcomputer  system  as  a 
demonstration  of  the  above-mentioned  possibilities.  The 
hardware  configuration  used  was  an  INTELLEC  8  microcomputer 
system. 


8 


The  architecture  of  the  microprocessor,  the  INTEL  8080, 
around  which  the  microcomputer  system  INTELLEC  8  is  built, 
consists  of  an  ALU  (arithmetic  logic  unit)  with  an  8  bit 
parallel  operation  and  with  the  option  of  either  binary 
coded  decimal  or  binary  arithmetic.  The  8080  has  a  16  bit 
address  capability,  which  allows  it  to  communicate  to  64K 
bytes  of  memory.  To  achieve  all  the  address  and  parallel 
data  capability,  the  8080  is  packaged  in  a  40  pin  package. 

In  developing  the  controller,  a  mathematical  model  for 
the  simulation  and  the  controller  is  developed  first.  These 
models  are  implemented  in  the  digital  simulation  language 
DSL/360.  The  control  part  then  is  coded  in  the  micro¬ 
computer  language  PLM  to  investigate  the  feasibility  of 
implementing  the  controller. 
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II.  MATHEMATICAL  MODEL 


There  are  six  equations  of  motion  of  a  submarine;  one 
for  each  degree  of  freedom.  Three  equations  are  for  linear 
motion  along  each  of  the  three  axes,  and  three  equations  are 
for  rotary  motion  about  each  of  the  three  axes  [Figure  ] ] . 

The  submarine  is  controlled  in  all  six  degrees  of  freedom 
by  a  set  of  movable  plane  control  surfaces,  movement  of 
ballast  water  and  the  propulsion  system.  The  set  of  movable 
plane  control  surfaces  consist  of  the  fairwater  planes,  the 
stern  planes  and  a  rudder. 

The  controller  for  depth  and  pitch  of  a  submarine 
described  in  the  following  sections  minimizes  the  sum  of  the 
errors  in  depth  and  in  pitch  together  with  the  control  effort. 
The  controller  counteracts  deviations  from  ordered  depth  and 
ordered  pitch  caused  by  external  forces  by  the  movements  of 
fairwater  and  stern  planes. 

The  simulated  submarine  is  a  hypothetical  one.  The 
dimensions  and  hydrodynamic  coefficients  of  this  submarine 
are  taken  from  Ref.  1.  The  implemented  controller  is  unique 
to  this  submarine,  but  the  underlying  idea  is  usable  in 
general.  The  basic  model  for  the  simulation  and  the 
controller  are  taken  from  Ref.  6. 

A.  EQUATIONS  OF  MOTION  FOR  A  SUBMARINE 

The  equations  of  motion  are  found  in  Ref.  2.  With  the 
assumption  that  the  principal  axes  of  inertia  of  the 
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Figure  1  -  Sketch  Showing  Positive  Directions  of  Axes,  Angles,  Velocities,  Forces,  and  Mo 


submarine  coincide  with  the  origin  placed  at  the  center  of 
gravity,  a  modified  form  of  the  equations  of  motion  is 
obtained  [Ref.  1] . 

In  Appendix  A  these  equations  are  linked  together  with 
the  nomenclature.  Positive  directions  of  translations  and 
rotations  are  in  the  directions  of  the  arrows  in  Figure  1. 

The  equations  of  motion  are  basically  of  two  forms 
[Ref.  4].  For  linear  motion: 

f  =  ma 

where 

f  =  force 
m  =  mass 

a  =  acceleration. 

For  rotary  motions: 

M  =  10 

where 

M  =  torque 

I  =  angular  moment  of  inertia 
0  =  angular  acceleration. 

The  actual  equations  are  of  the  form: 

£  (ma)  =  J  f  and  £  (10)  =  J  M  . 

The  equations  contain  forcing  terms  for  the  planes,  the 
rudder,  the  propeller  and  the  trimming  of  ballastwater .  The 
input  for  the  propeller  is  a  function  obtained  by  curve 
fitting.  There  are  different  coefficients  for  different 
propeller  modes  (forward,  backing,  acceleration,  slowing 
down) .  Only  the  coefficients  for  forward  motion  at  constant 
speed  were  used  in  this  investigation. 

The  manipulation  of  these  equations  to  get  the  expressions 
for  u,  v,  w,  p,  q,  r  are  shown  in  Appendix  A. 
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B.  DEPTH  AND  PITCH  CONTROL  USING  OPTIMAL  CONTROL  THEORY 
To  use  the  equations  of  motion  in  the  controller,  the 
nonlinear  equations  were  linearized.  The  linearization  and 
the  assumptions  made  are  shown  in  Appendix  A.  With  the 
resulting  equations  involving  w,  the  vertical  velocity 
component,  and  qr  the  angular  velocity  about  the  body  y-axis, 
the  problem  is  reduced  to  two  degrees  of  freedom.  The 
controls  and  6^  correspond  to  deflections  of  the  stern 
plane  and  bow  plane,  respectively. 


m'-z: 

w 

-*.z: 

q 

w 

'vy* 

(m'+Z  )u 
q  3 

w 

-M \/l 
w 

I'-M! 

y  qj 

q 

vs/*2 

u  M'A 
a  q 

q 

+ 


uJz’SsA  u2aZ'&b/i 

u2M'4b/^ 


More  generally,  the  linearized  equations  written  in  state 
equation  form  are 


(2)  x  =  Gx  +  Bu  , 


where 


13 


X 


G 


B 


0  10  0 

0  g22  0  g24 

0  0  0  1 

0  g42  0  g44 

0  0 

b21  b22 

u  = 

0  0 

b41  b42 


with 

g22  =  tCi;-M*)Zw  +  M^Z^]ua/(Z-DET) 

g24  =  (ra'  +Zq>  +  M^]ua/DET 

g42  =  +  -Z£)]ua/U2-DET) 

g. .  =  [M! (m’  +  Z’ )  +  M'  (m'  -  Z! ]u  /U-DET) 
r44  w  q  q  w  a 

b21  =  t(X;-MVZis  +  M^Z(-]u2/(£-DET) 
b22  =  [(Iy"M4)Z^b  +  M^]u2/(i’DET) 
b41  =  +  (m'  -Z^)M^]u2/(Z-DET) 

b42  =  [M*Z5b  +  (m1  -Z*,)M^]u2/Z2-DET) 


DET  =  (m1  -Z^)  (ly  -  M^)  -  Z!H! 
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Referring  to  Ref.  3,  the  linear  tracking  model  was  used 
to  solve  the  problem  because  the  desired  value  of  the 
state  vector  is  not  the  origin.  The  reference  vector  in 
this  case  is 


with  r^  =  ordered  depth  and 
r^  =  ordered  pitch, 
both  are  considered  to  be  constant. 
Let 

y  -  x  -  r 

(3)  y  =  x  -  r 

and  substituting  (2)  into  (3)  yields 


(4) 


y  =  Gx  +  Bu  -  r 


Realizing  that 


r 


0 


and 


Gr  = 


0 

0 

0 

0 
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equation  (4)  can  be  written 


y  =  G[x  -  r]  +  Bu 


(5) 


y  =  Gy  +  Bu 


The  performance  measure  to  be  minimized  may  be  expressed  as 


(6) 


’o 


where  the  final  time  t^  is  fixed  and  y(t^)  is  free.  Q  is  a 
real  symmetric  positive  semi  definite  weighting  matrix  and 
S  is  a  real  symmetric  positive  definite  weighting  matrix. 
Then  (5)  and  (6)  define  a  linear  regulator  problem. 

The  Hamiltonian  is  given  by 


H (y (t) f u(t) >p (t) ,  t)  =  j  [yT(t)Qy(t)  +  u(t)TSu(t)] 


+  p ( t) TGy ( t)  +  p(t)TBu(t) 


The  costate  equation  is 


(7) 


P*(t) 


f|  =  -Qy(t)  -  GTp*(t)  .  1 


Optimality  is  indicated  by  asterisks. 
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The  algebraic  relation  which  has  to  be  satisfied  is 
0  =  |S  =  Su* (t)  +  BTp*(t) 
which  gives  the  control  vector 
(8)  u*(t)  =  -S-1BTp*(t) 

whose  elements  are  the  control  plane  deflections  g  and  ^ 
Substituting  (8)  in  the  state  equation  (5)  and  forming  the 
augmented  matrix  yields 


y*(t) 

i  -1  T 

G  1  -BS  Bx 

y*(t) 

P*(t) 

l 

O 

1 

i _ 

P*(t) 

Reference  3  gives  the  boundary  condition  as  p*(t^)  =  0. 
Defining 


(10) 

p*  (t) 

=  K(t)y* 

(t) 

[Ref.  3] 

(ID 

p*(t) 

=  K(t)y* 

(t) 

+  K  (t)  y*  (t) 

and  combining  (7)  and  (11)  gives 

-Qy*(t)  -  GTp*(t)  =  K(t)y* (t)  +  K(t) y*(t)  . 


17 


Solving  for  K(t) ,  substituting  y*(t)  from  (9)  and  replacing 
p*(t)  by  (10)  yields 

K(t) y* (t)  =  -K (t) Gy* (t)  +  K(t)BS_1BTK(t)y*(t) 

-  Qy* (t)  -  GTK(t)y*(t) 

which  must  be  satisfied  for  arbitrary  y* (t) .  Hence 

(12)  K(t)  =  -K(t)G  +  K(t)BS_1BTK(t)  -  Q  -  GTK(t)  . 

The  solution  to  equation  (12) ,  the  Riccati  equation,  is  the 
matrix  of  the  optimal  parameters  for  the  defined  mathematical 
model  based  on  the  quadratic  performance,  measure  (6) . 
Equations  (8)  and  (10)  give  the  final  expression  for  the 
control  vector  in  terms  of  the  matrix  K 

u*(t)  =  -S~1BTK(t)y* (t) 
or 

(13)  u*(t)  =  S-1FTy*(t) 
where 

F (t)  =  -K(t)B  . 

K ( t)  is  a  symmetric  4x4  matrix,  so  that  ten  differential 
equations  for  K(t)  must  be  solved  to  determine  the  optimal 
gains  for  u*(t) .  The  solution  process  is  described  in 
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section  III.A.1.  From  Ref.  3  K(tf)  =0.0.  The  matrices  Q 
and  S  are  weighting  matrices  and  have  the  form 


where  w^  is  the  weight  for  the  error  in  depth,  w2  is  the 
weight  for  the  error  in  pitch  and  c1  and  c2  are  the  weights 
for  the  control  input. 

Since  only  the  value  of  K(t^)  is  known,  equation  (12) 
has  to  be  solved  backwards.  By  multiplying  the  right  side 
of  (12)  by  -1.0  and  reversing  the  order  of  integration, 
K(tf)  becomes  the  initial  condition. 
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III.  COMPUTER  IMPLEMENTATIONS  OF  THE  MATHEMATICAL  MODEL 


The  programming  language  for  the  simulation  and  the 
controller  of  the  submarine  was  DSL/360,  a  digital  simulation 
language.  The  programs  were  run  on  the  IBM  360/67  computer 
at  the  Naval  Postgraduate  School. 

For  the  implementation  of  the  controller  on  an  INTELLEC  8 
microcomputer,  the  programming  language  PLM  was  used.  The 
PLM  programs  were  compiled  with  a  cross-compiler  on  the 
IBM  360/67. 

A.  DSL/360  PROGRAMS 

The  DSL  programming  part  was  divided  into  two  programs. 
The  first  program  (GAIN)  performs  the  computations  which 
have  to  be  done  only  once.  The  results  were  transferred 
to  a  datacell  for  later  use  by  the  second  DSL/360  program 
(SIMUL) ,  the  actual  simulation  of  the  submarine.  Figure  2 
shows  a  general  flow  diagram  for  both  DSL/360  programs. 

Figure  3  gives  a  graphical  representation  of  the  different 
timing  intervals.  A  communications  interval  in  DSL/360 
is  an  interval  after  which  information  is  returned  to  the 
user.  In  the  case  of  these  two  programs  the  sample  region 
is  called  in  communication  time  intervals  to  collect  data 
points  for  the  graphical  output.  The  length  of  the 
communication  interval  is  determined  by  the  program  control 
variable  DELS. 
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CONNECTOR  FOR  ENTRY 
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Figure  2 


21 


Figure  2  (Continued) 
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A  simulation  interval  is  the  unit  of  time  for  a  partic¬ 
ular  integration  routine  to  complete  an  integration  step. 

The  length  of  the  simulation  interval  is  determined  by  the 
program  control  variable  DELT. 

A  calculation  interval  is  a  sub interval  of  the  simula¬ 
tion  interval  whose  length  is  determined  by  the  integration 
method  used.  In  these  two  DSL/360  programs  the  variable 
length  step  Runge  Kutta  (RKS)  integration  routine  was  used 
with  a  calculation  interval  of  0.25*DELT.  Reference  5  gives 
a  more  detailed  description  of  DSL/360. 

1.  DSL/360  Program  GAIN 

The  program  GAIN  uses  the  library  routine  MINV  to 
get  the  inverse  of  the  matrix  A  to  solve  the  equations  of 
notation  Cl) » 


A«  DOT  =  Z 

by  the  program  SIMUL. 

The  kernel  of  the  program  GAIN  is  the  equation  (12) 

K(t)  =  — K (t) • G  +  K(t)BS_1BTK(t)  -  Q  -  GTK(t)  . 

The  matrices  G,  B,  R  and  Q  are  constant.  Under  the 
assumption  that  the  plant  is  completely  controllable, 

K ( t)  K  as  tf -►<»  [Ref.  3].  But  even  if  the  control  interval 
is  finite  and  long  with  respect  to  the  decay  time  of  the 
optimal  gains,  it  is  possible  to  use  the  fixed  control  law 
(see  Figs.  4.1  to  5.8). 


24 


To  compute  the  steady  state  gains  K,  either  the  Riccati 
equation  (12)  is  integrated  or,  setting  K(t)  =  0,  the 
algebraic  equation 

0  =  -KG  -  GTK  -  Q  +  KBS-1BTK 

is  solved.  The  method  of  integrating  equation  (12)  was 
chosen  because  of  the  graphical  output  (Figs.  4.1  to  5.8) 
which  shows  the  behavior  of  K(t)  for  t The  gain  is  also 
a  function  of  the  submarined  speed  ua*  Therefore  different 
gains  were  computed  for  a  representative  set  of  speeds  u&. 

This  approach  may  result  in  suboptimal  gains  when  speed  u_ 

3. 

deviates,  e.g. ,  in  a  turn  or  during  a  change  of  depth  (see 
Figs.  9.5  and  10.5).  But  even  in  the  maneuver  of  a  turn 
(Fig.  10.5)  the  temporary  drop  in  speed  ua  occurs  over  a  time 
period  which  is  long  compared  to  the  decay  time  of  the  optimal 
gains . 

Together  with  the  computations  for  different  speeds  u&, 
different  weights  w^  and  w ^  were  chosen  to  have  the  possi¬ 
bility  of  adjusting  to  environmental  effects  (seastate,  depth 
of  water,  on  periscope  depth,  speed,  etc.).  The  gains  for 
these  different  combinations  were  computed  in  repeated 
DSL/360  runs. 

The  values  chosen  for  this  hypothetical  submarine  were: 


i) 

Speed  ua 

from  1.689  ft. /sec. 

to  30.402 

ft. /sec.  in  increments 

of  1.689 

ft. /sec.  =  l.Okn 

ii) 

Weight  w^ 

:  0.8,  0.9,  1.0,  1.1,  1.2 

iii) 

Weight 

:  2400.0,  2700.0,  3000,0,  3300.0, 

3600.0 
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These  values  depend  on  the  specific  type  of  submarine  and 
have  to  be  determined  experimentally  or  by  a  simulation  for 
that  specific  submarine. 

The  steady  state  gains  were  assumed  after  70  sec. 
for  small  ua's  and  sooner  for  larger  ua's  (see  Figs.  4.1 
to  5.8).  The  length  of  each  run  was  set  to  200  sec.  The 
sets  of  gain  values  were  collected  in  a  five-dimensional 
array  FF  which  was  printed  and  transferred  to  a  datacell 
after  the  last  run. 

2.  DSL/360  -  Program  SIMUL 

Program  SIMUL  consists  of  the  main  program  and  the 
subroutines  CNTRL,  PLGEN,  and  PRIME.  The  program  is 
centered  around  the  equations 

(1)  DOT  =  A_1-Z 

and 

(13)  U(t)  =  S-1FTyMt)  . 

To  use  these  equations  the  inverse  of  A  and  the  matrix  FF 
containing  all  F(I,J)  values  have  to  be  read  from  the 
datacell  at  the  beginning  of  the  program. 

With  the  integrated  elements  of  DOT  the  values  of 
the  variables  PITCH  and  DEPTH  are  computed.  These  variables 
are  controlled  by  the  automatic  controller.  To  manipulate 
these  variables,  the  movements  of  fairwater  and  stern  planes 
have  to  be  simulated.  For  turns  about  the  vertical  axes  the 
rudder  is  simulated,  too. 
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F ( 1 . 1 )  US  TIME 
un  =  5»0?  FT/SEC 


c* 


VSCRLl."C.  40  UNITS/INCH  PLOT  N0«  L 

FIGURE  4.1 
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F ( 1 ..  2 )  US  TIME 
UR  =  5»0?  FT/SEC 


'^SCflLF -40.  00  UN  ITS/ INCH  RUN  MO.  I 

VSCRLE-C.40  UNITS/ INCH  PLOT  NO. 2 


FIGURE  4.2 
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FC  2, 1)  US  TIME 
UR  =  5*  0?  FT/SEC 

CJ 

c> 


VSCRLE>5. 00  UN ITS/ INCH  PLOT  NO. 3 

FIGURE  4.3 
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F(  2,  2!  US  TIME 
UR  =  5-0?  FT/5EC 


;kio  ; 

K5CRLE--  40.00  UHIT5/INCH  RUN  HO.  L 

V5CRLF.-3*  OO  UNITS/INCH  PLOT  NO,  4 


Figure  4.4 
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F (  3 .•  1 )  US  TIME 
UR  -  5=0?  FT/SEC 


'  KiO  / 

^5CRLur-40*  l Q  UNITS'INCH  RUN  NQ.  L 

y5CRLE-20.0C  UNITS'INCH  PLOT  NQ,5 

Figure  4 . 5 


time 
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F { 3 . 2 )  U3  TIME 
UR  =  5.0?  FT/SEC 


X5CRLE-40.00 

■m  i  1 

UNITS? INCH 

RUN  NO-  i 

V5CRL5-2Q.00 

UN  ITS/ INCH 

PLOT  NO*  S 

Figure  4 . 6 
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Ft  4. 1)  US  TIME 
UR  =  5-0?  FT/SEC 


b 


XSCRLE-40.  00  UNITS/' INCH 

''SCRLE--.-40C.00  UNITS/IHCH 


RUN  NO.  L 

plot  no. ? 


FIGURE  4.7 


33 


F( 4 >  2 J  US  TIME 
UR  -  5-  0?  FT/SEC 


'''5r-RLE“40u»  00  LlNITS/IHCh  PLOT  NO.  5 

FIGURE  4.8 
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F( 1, 1)  US  TIME 
UR  -  25. 35  FT/SEC 


XSCRLE--40. 00  UNITS/INCH  RUN  NQ.  i 

VSCRLE=0. 40  UNITS^INCH  PLOT  NO. 1 


FIGURE  5.1 


TIME 


35 


F(1,2J  US  TIME 
UR  =  25.35  FT/SEC 


;kio  i  1 

XSCflLE--=40.00  UNITS/INCH  RUN  NO.  1 

VSCRLE-0.40  UNITS/INCH  PLOT  N0>2 

FIGURE  5.2 


TIME 
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F(2,l)  US  TIME 
UR  -  25.35  FT/SEC 


i) 


~~i  i 

■3.33  i.OC 

KSCRLF.----40.00 
ySCRLE-l. 00 


1 - i - 

3.QQ  12-03 

ikio  j  1 

UHITS/INCH 

UNITS/INCH 

FIGURE  5.3 


H - 1 

Is).  00  13.33 

RUN  NG*  i 
PLOT  NO. 3 


TIME 
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F( 2. 2)  U5  TIME 
UR  =  25.35  FT/SEC 


'K1Q  l  1 

XSCnLE-40.00  UHITS/IHCH  RUN  NO, i 

ySCRLE-2.00  UNITS/INCH  PLOT  NO, 4 

FIGURE  5.4 


TIME 
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FC3.11  US  TIME 
UR  =  25.35  FT/SEC 


& 


v  !K10  I  1 

KSCRLE--40. 00  UNITS/INCH  RUN  HO.  1 

V3CRLE-20. 00  UNITS/INCH  PLOT  HQ. 5 

FIGURE  5.5 


TIME 


39 


F (  3 •  2 )  US  TIME 
UR  =  25-35  FT/SEC 


4-3: 


s.oo 


’.KIO  } 


12. 20 

1 


IS.  00 


20.33 


TIME 


KSCRLE-40. 00 
¥SCRLE>20. 00 


UNITS/ INCH 
UNITS/INCH 

FIGURE  5.6 


RUN  NQ-1 
PLOT  NO. 6 
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F ( 4,1)  US  TIME 
UR  =  25. 35  FT/SEC 


*SCflLE----40.00  UNITS/INCH  RUN  HO.  L 

VSCRLE--80. 00  UNITS/INCH  PLOT  NO.? 


FIGURE  5.7 
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F(4,2)  US  TIME 
UR  =  25-35  FT/SEC 


'-‘SCPLr  -50. 00  UNITS'INCH  PLOT  NO. 6 

FIGURE  5.8 


TIME 
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To  change  the  various  command  inputs,  data  cards 
with  the  given  values  are  read.  The  first  variable  read  (TT) 
is  the  time  of  the  next  change  of  command  followed  by  the 
numerical  values  of  the  command  variables.  The  format  is 
shown  in  Fig.  6.  At  least  one  data  card  has  to  be  used. 

If  no  changes  in  command  variable  values  are  intended,  the 
specified  time  TT  on  the  first  card  has  to  be  one  not  in 
the  range  between  0.0  and  FINTIME,  the  end  of  the  simulation. 

At  the  end  of  section  III. A. 2. c  some  graphs  of 
simulation  runs  are  shown. 

a.  Subroutine  CNTRL 

This  subroutine  represents  the  automatic 
controller  and  is  called  with  the  following  arguments: 

IW1,  the  choice  of  weight  for  the  error  in  depth; 

IW2,  the  choice  of  weight  for  the  error  in  pitch? 

Cl  ,  the  weight  for  the  control  input  D(l); 

C2  ,  the  weight  for  the  control  input  D(2); 

ZORD,  the  ordered  depth? 

PORD,  the  ordered  pitch. 

Beside  these  command  inputs,  the  measured  inputs  DEPTH,  PITCH 
and  speed  ua  are  passed  through  the  argument  list. 

The  speed  ua  is  truncated  to  an  integer  IU  (in 
increments  of  1.689  ft. /sec.).  With  IU,  IW1  and  IW2  as 
subscripts  for  the  array  FF  the  gain  values  F(I,J)  are 
extracted.  Together  with  the  errors  between  actual  depth 
and  ordered  depth,  actual  pitch  and  ordered  pitch,  rate  of 
change  in depth  and  pitch  returned  from  routine  PRIME,  the 
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f F10.4 

F10.4 
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31  40 

41  50 
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71  75 

76  80 

TT 
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ZORD 

FORD 

DRORD 

Cl 

C2 

IW1 

IW2 

TT 

UC 

ZORD 

PORD 

DRORD 

Cl 

C2 

IW1 

IW2 


Time  of  next  change  of  command  variables 
Command  speed 
Ordered  depth 
Ordered,  pitch 

Ordered  deflection  of  rudder 
(for  the  simulation  of  turns) 

Weight  of  control  input  for  sterm  plane 
(range  8.0  to  12.0) 

Weight  of  control  input  for  fairwater  plane 
(range  8.0  to  12.0) 

Choice  of  weight  for  error  in  depth 
(range  1,  2 ,  . . . ,  5) 

Choice  of  weight  for  error  in  pitch 
(range  1,  2,  . ..,  5) 


FIGURE  6 


A  A 


control  outputs  D(l)  and  D ( 2 )  are  computed  with  equation  (13) 
These  outputs  are  passed  back  to  the  simulation  over  the 
COMMON- area  Rl. 

Graphical  outputs  for  a  change  in  speed  (Figs. 

8.1  to  8.5),  a  turn  (Figs.  9.1  to  9.5)  and  a  change  in 
depth  (Figs.  10.1  to  10.5)  are  shown  at  the  end  of  Section 
III. A. 2.  All  runs  used  error  channel  limiter  and  plane 
deflection  limits.  Each  set  of  graphical  outputs  consists  of 
a  graph  for 

1)  stern  plane  angle  vs.  time 

2)  fairwater  plane  angle  vs.  time 

3)  depth  vs.  time 

4)  pitch  vs.  time 

5)  speed  ua  vs.  time. 

For  a  60  ft.  change  in  depth  without  constraints, 
the  results  are  shown  in  Figs.  11.1  to  11.5.  The  large 
deflections  of  the  planes  are  not  desirable  and  in  practice 
are  also  not  achievable.  The  reason  for  this  behavior  is 
a  non-optimal  distribution  of  the  weights  for  the  errors  and 
for  the  controls.  The  magnitude  of  the  control  outputs  has 
to  be  smaller.  Two  possibilities  exist  to  achieve  this 
result:  either  redistribution  of  the  weights  in  the 

weighting  matrices  or  limitation  of  the  control  inputs.  In 
this  case  the  control  inputs  "error  in  depth"  and  "error  in 
pitch"  have  been  limited  to  4.0  ft.  and  10.0  degrees. 

The  latter  approach  is  a  compromise  between  the 
decrease  of  the  control  outputs  and  a  good  response  for  small 
errors  in  depth  and  pitch. 
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b.  Subroutine  PRIME 

Instead  of  taking  the  rate  of  range  in  depth 
and  in  pitch  from  the  simulation  (ZODOT  and  PIDOT) ,  this 
routine  is  used  to  keep  the  controller  as  independent  as 
possible  from  the  simulation. 

To  get  the  rate  of  change,  a  linear  combination 
of  discrete  orthogonal  Legendre  polynomials  O^fx)  has  been 
used  to  fit  a  polynomial  curve  in  the  least  square  sense  to 
a  set  of  points  with  evenly  spaced  abscissas  [Ref.  7].  The 
discrete  Legendre  polynomials  for  a  second  degree  polynomial 
are 


°0n  1 


°m  =  1  -  2<k> 


°2n  =  1  -  6(K>  +  • 


Using  the  properties  of  orthogonal  polynomials 


JQ  °in(x)  °jn(x)  =  0  '  1  *  3 


the  normal  equations  for  the  least  squares  fit  are 


(21) 


Jo0-'*’  “ 


I.  °1»  (,)  0 

x=0 


o  y  o*{x) 

x=0  2n  . 


Jo  °0n(x)y* 

1  °ln(x)yx 


x=0 


l-  02n(x)^x 


Lx=0 
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2 

where  the  yx  represent  the  sampled  data  points.  The 
linear  combination  of  Legendre  polynomials  which  gives  the 
least  squares  fit  is 


(22) 


F(x) 


=  a0°0n(x)  +  al°ln(x)  +  a2°2n(x) 

=  a0[l]  +  a,  1 1  -  2  f]  +  a2[l-6  |  + 


6 


x(x  -  1) 

nTn  -1)  J 


Transforming  the  set  of  evenly  spaced  points  x"  into  integer 


(23)  x  =  -  -h~  —  with  h  =  xi+i“x|  =  constant 

h  =  DELS  in  program  implementation. 


Taking  five  points  (n=4)  and  evaluating  the  polynomial  at 
the  third  point  (Z=0.4)  the  values  of  the  orthogonal 
polynomials  are 


and 

1  00.2(X)  =  5  l  0142.(x)  =  2.5  l  0242(x)  =  3.5 

x=0  x=0  x=0 


2 

X  denotes  an  integer  0 ,  1 ,  2 . 
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From  (21)  then,  using  the  above  numerical  values 

ao  =  I  <y0  +  +  y2  +  y3  +  y4) 

al  =  k  (2y0  +  yl  -  y3  -  2y4> 
a2  "  7  (2y0  "  yl  “  2y2  -  y3  +  2y4> 
and  (22)  becomes 

F (x)  =  aQ  +  ax  +  a2  -  (-^  +  2a2)x  +  -y-  x2 
The  derivative  of  F(x)  is 

F* (x)  =  -  j  ai  ”  2a2  +  a2x  • 


Replacing  a^,  a^,  a2  and  using  (23)  gives 

F'(x')  =  -  ^  (54y0  -  13y1  -  40y2  -  27y3  +  26y4) 

+  JK  (2y0  -  yl  "  2y2  -  y3  +  2Y4>  <*'  +2-°) 

If  x*  is  the  time  of  sampling,  h  the  sampling  interval  and 
yx  the  sampled  data  points,  then  F'(x')  is  the  rate  of 
change  in  F  at  x* . 

c.  Subroutine  PLGEN 

This  routine  is  called  three  times  in  the 
derivative  region  for  the  simulation  of  the  movements  of 
fairwater  and  stern  planes  and  the  rudder.  The  time  between 
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the  repeated  calls  is  determined  by  the  length  of  the 
calculation  interval.  The  routine  is  called  with  the 
following  arguments: 

(1)  deflection  of  the  plane  or  rudder; 

(2)  difference  between  deflection  and  ordered 
deflection  from  the^ controller ; 

(3)  increment  or  decrement  of  deflection  computed 
in  the  previous  call; 

(4)  rate  of  change  in  the  plane  movement 7 

(5)  length  of  time  interval,  in  this  case  DELT  *  0.25, 
the  length  of  the  calculation  interval. 

The  increment  in  the  plane  or  rudder  deflection  computed 
in  the  previous  call  is  added  to  the  current  plane  or  rudder 
deflection  unless  the  increment  is  smaller  than  the  dead 
zone  of  t  0.5  degrees.  A  limiter  limits  the  maximum  plane 
or  rudder  deflection  to  36.0  degrees.  A  flow  diagram  is 
shown  in  Figure  7. 


49 


Figure  7 


50 


Figures  8.1  to  8.5 


uc  =  25.34  ft. /sec. 


change  of  uc  from  25.34  ft. /sec.  to  15.21  ft. /sec. 
at  TIME  =  50.0  sec. 
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STERN  PLANE  ANGLE  US  TINE 


Nil  UN  I T  =  1  SEC.  Y : 1  UNIT=1  OEGR 


XSCPlLE=eO. DO  UNITS'IhCH  RUN  HQ.  1 

V5CRLE--12. 00  UNITS/INCH  PLOT  NO.l 


FIGURE  8.1 
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FRIRRRTER  PLRNE  RNGLE  US  TIME 
Xil  UNIT=1  SEC,  V : 1  UNIT=1  DEGR 


IK10  J  1 


TIME 


XSCRLE=80.00  UNITS/ INCH  RUN  NO. 1 

v5CRLE~i2. 00  UNITS/INCH  PLOT  NO. 2 


FIGURE  8.2 
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DEPTH  US  TIME 

X:.t  UNIT=1  SEC.  V:  L  UNIT-1  FT 


“SC-ALE'l.OO  UNITS'INCH  PLOT  HO.  3 


FIGURE  8.3 


TIME 


54 


PITCH  US  TIME 

X:1  UNIT=1  SEC,  V:1  UNIT=1  DEGR 


KSCRLE-80.00  UNITS/INCH  RUN  HO. 1 

VSCRLE-2. 00  UNITS/INCH  PLOT  NO. 4 


FIGURE  8.4 
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SPEED  UR  US  TIME 

X:1  UNI T=1  SEC,  V:1  UNIT=1  FT/SEC 


KSCRLE-80.00  UNITS/INCH  RUN  N3-1 

YSCRLE=2- 00  UNITS/INCH  PLOT  NO. 5 


FIGURE  8.5 


^TTME 
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Figures  9,1  to  9.5 


uc  =  25.34  ft. /sec. 


turn  with  a  rudder  deflection  of  20.0  degrees 
at  TIME  =  50.0  sec. 


57 


STERN  PLRNE  RNGLE  US  TIME 
X : 1  UNIT=1  SEC,  V:i  UNIT=1  DEGR 


YSCRLE-L2. 00  UNITS/ INCH  PLOT  NO.l 

FIGURE  9.1 
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FRIRWRTER  PLRNE  RNGLE  US  TIME 
X:1  UN  I T  =  1  SEC,  V: 1  UNIT=1  DEGR 


;kio  )  1 

XSCRLE=80.00  UNITS/INCH  RUN  HQ. 1 

Y5GRLE-12.00  UNITS/INCH  PLQT  HQ. 2 

FIGURE  9.2 
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DEPTH  US  TIME 

X : 1  UN  I T  =  1  SEC,  V: 1  UN  I T  =  1  FT 


;kio  i  1 

XSCRLE=80.00  UHIT5/INCH  RUN  HQ. 1 

VSCRLE^l.OO  UHITS/1NCH  PLOT  HQ. 3 

FIGURE  9.3 
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PITCH  US  TIME 

X:1  UNI T=i  SEC,  Y:i  UNIT=1  OEGR 


:xio  )  L 

XSCRLE=80. 00  UNITS/INCH  RUN  NcM 

VSCRLE=2. 00  UNITS/INCH  PLCT  NC*4 

FIGURE  9.4 
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SPEED  UR  US  TIME 

X:1  UN  I T =1  SEC,  V: 1  UNIT=1  FT/SEC 


PSCRLE--2. 00  UNITS/ INCH  PLOT  NO. 5 

FIGURE  9.5 
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Figures  10.1  to  10.5 


uc  =  25.34  ft. /sec. 


change  of  depth  by  60  ft.  at  TIME  =  100.0  sec. 

with  error  channel  limiter  and  plane  deflection  limits 
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STERN  PLRNE  RNGLE  US  TIME 
X:1  UNIT=1  SEC,  V: 1  UNITU  DEGR 


:sio  i  1 

XSCRLE=80. 00  UNITS/INCH  RUN  HO. 1 

V5CRLE--12. 00  UNITS/INCH  PLOT  NO.  i 

FIGURE  10.1 
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FAIRWATER  PLANE  ANGLE  US  TINE 
X: L  UNIT=1  SEC,  Y:1  UNIT=1  DEGR 


XSCRLE=;60.00  UNITS'INCH  RUN  HQ. 1 

V3CRLE-12.00  UNITS/INCH  PLOT  NQ.2 


FIGURE  10.2 
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DEPTH  US  TIME 

X21  uni t=i  SEC?  ya  unit=i  ft 


X5CRLE-80. 00  UNITS/INCH  RUN  H0.1 

V5CRLE-12.00  UNITS/INCH  PLOT  NO. 3 


FIGURE  10.3 
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PITCH  US  TIME 

X:1  UNI T=1  SEC,  V: 1  UNIT=1  DEGR 


;kio  i  1 

*SGflLE=80.0G  UNITS/INCH  RUN  NO, I 

VSCRLE--2. 00  UNITS/INCH  PLOT  N0.4 

FIGURE  10.4 
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SPEED  UR  US  TIME 

X: 1  UNIT=1  SEC,  V:1  UNIT-1  FT/SEC 


ySCRLE~2.  00  UNITS/INCH  PLOT  NO.  5 

FIGURE  10.5 


TIME 
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Figures  11.1  to  11.5 


uc  =  25.34  ft. /sec. 


change  of  depth  by  60  ft.  at  TIME  =  100.0  sec. 
without  error  channel  limiter  and  without 
plane  deflection  limits 
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STERN  DLRNE  ANGLE  US  TIME 
X:.L  UNIT-1  SEC,  V:.l  UNIT*!  DEGR 


VSCRLF-30.00  JN ITS/ INCH  PLOT  HC. I 

FIGURE  11.1 


TIME 
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FRIRMRTER  DLRNE  RNGLE  US  TIME 
X:  1  UNIT=1  SEC.  V: 1  UN  I T=1  DEGR 


iWM  i  1 

'KSCALE'30*  00  UNITS/INCH  RUN  NO*  L 

‘J5CPLE -100*  Q0  UNITS/INCH  PLOT  NO.  2 

FIGURE  11.2 
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DEPTH  U5  TIME 

Ms  1  UNI T=1  5EG.  V: 1  UN  I  T=1  FT 


:KiO  . 

K5CRLE--30.00  UNITS/INCH  RUN  HO.  1 

VSCRLF--  20.00  UNITS/INCH  PLOT  NO- 3 


FIGURE  11.3 


TIME 
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CO  CO 


PITCH  US  TIME 

K:i  UNIT-1  SEC-  Y: t  UNI T=1  DEGR 


:klo  j 

CRLF-SG*  00  UHITS/IHCH  RUN  HO. 1 

CPLC-3*  00  UHITS/IHCH  PLOT  HO.  4 


FIGURE  11.4 


73 


SPEED  Ufi  US  TIME 

Mil  UNIT=1  SEC.  V : 1  UN  I  T=1  FT/SEC 


\klo  ; 

:K5CRLF---30.00  UNITS/INCH  RUN  HO.  1 

YSCRLE---4.00  UN  IT  5/ INCH  PLOT  NO.  5 


FIGURE  11.5 
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B.  PLM  -  PROGRAM 


1 .  Data  Formats  and  Conversions 

The  hardware  of  the  INTEL  8080  microcomputer  is 
restricted  to  integer  arithmetic.  The  basic  wordlength 
is  one  byte.  Floating  point  arithmetic  is  possible  by  means 
of  software  routines.  The  control  algorithm  uses  routines 
for  multiplication  and  addition  in  floating  point  format. 

The  floating  point  representation  has  the  following  form. 

MANTISSA  EXPONENT 

I - 1 - 1 

0  15  16  23 

The  mantissa  is  normalized,  the  high  order  bit  in  the 
exponent  represents  the  sign  of  the  number.  The  lower  half 
of  the  remaining  range  of  the  seven  bits  (0  to  2^  -  1)  is 
used  for  exponents  less  than  0,  the  upper  half  (26  to  27  -  1) 
is  used  for  exponents  equal  or  greater  than  0.  The  smallest 
magnitude  which  can  be  represented  in  floating  point  format 
is  1.8446  •  10  and  the  largest  representable  number  is 
9.2233  •  1018. 

The  format  of  the  input  and  output  of  the  control 
algorithm  is  contained  in  a  16  bit  word  which  represents 
an  integer  in  the  range  from  0  to  65535.  The  range  between 
0  and  32768  represents  negative  values,  the  range  over  32768 
represents  positive  values.  The  contents  of  this  16  bit 
word  is  assumed  to  come  from  an  A/D-converter  and  to  go  to 
a  D/A-converter .  The  range  of  the  16-bit  integer  representa¬ 
tion  corresponds  to: 
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depth: 


0  to  32767  inches  =  32768  to  65535 

(0  to  32767  not  used) 
pitch:  16.0  degrees  to  +16.0  degrees  =  0  to  65535 

1.0  degrees  =  2048 

4.883  •  10  ^  degrees  =  1 

plane  deflection:  -36. GO degrees  to  +36.0  degrees 

=  0  to  65535 

1.0  degrees  =  910  approx. 

1.0987  degrees  =  1 

There  are  two  routines  for  the  necessary  conversions: 

•  routine  FLOAT  converts  from  integer  to  floating  point 
format. 

•  routine  INTEG  converts  from  floating  point  to  integer 
format. 

The  use  of  the  time  consuming  arithmetic  software  routines 
has  been  limited  to'  the  minimum  necessary.  All  computations 
which  could  be  done  outside  the  PLM-program  are  done 
beforehand. 

The  results  then  are  used  in  the  PLM-program.  Where 
possible,  comparisons  and  arithmetic  manipulations  have  been 
done  with  the  faster  integer  format.  This  results  in  an 
execution  time  of  ca.  0.12  sec.  for  one  pass  through  the 
control  algorithm.  It  is  important  that  the  time  for  one 
pass  does  not  exceed  0.2  sec.  Otherwise  the  sampling  interval 
for  depth  and  pitch  becomes  too  large  and  the  controller 
does  not  work  properly.  This  observation  results  from  the 
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DSL/360  simulation.  The  PLM-version  of  the  controller 
consists  only  of  one  pass  through  the  algorithm  and  only 
of  two  sets  of  gain  values  in  the  F  array.  This  was 
considered  to  be  sufficient  for  demonstrating  the  feasibility. 
In  an  actual  implementation  the  F-array  had  to  have  the  same 
number  of  elements  as  in  the  DSL-programs .  Furthermore,  the 
controller  would  have  to  work  in  an  infinite  loop  with  fixed 
intervals  for  reading  the  input  values  from  the  input  ports. 
This  sampling  would  replace  parts  of  the  input  routines, 
conversion  routines  and  input  statements. 

2 .  Error  Analysis: 

A  floating  point  mantissa  consists  of  16  binary 
digits,  or  approximately  five  decimal  digits.  Starting  with 
the  assumption  that  the  magnitude  of  the  error  connected  with 
each  number  consists  only  of  roundoff  errors,  then  the 
maximum  error  per  number  can  be  j  •  2  =2 

As  an  example,  the  absolute  error  bounds  are  computed 
for  D ( 1)  and  Y(2)  (see  Appendix  B) .  The  results  are: 

|  De  (1)  |  £  |  C6  •  Y  ( 3)  •  F(4,l)  (A6  +  61  +  al  +  a2  +  a3  +  ml  +m2) 

+  F  ( 3 , 1)  •  Y(4)  (f5 +Ye(2)  +  m3  +  al  +  a2  +  a3) 

+  F  (2 , 1)  •  Y(2)  (f3  +Ye(4)  +  m4  +  a2  +  a3) 

+  F  ( 1 , 1)  •  Y(l)  (fl  +  62  +m5  +  a3)  | 


and 


77 


|  Ye  (2)  |  <  Cl  •  dl(Al  +  51  +  ml  +al  +  a2  +a3  +a) 
+  C2  •  d2  (A2  +  62  +  m2  +  al  +  a2  +  a3  +  a) 
+  C4-d4(A4+64+m4+a2+a3+a) 

+  C5*d5(A5+55+m5+a3+a) 

+  C3  •  d3  (A3  +  53  +  m3  +a) 


where 

Ai  =  roundoff  error  in  coefficient  for  Legendre  Polynomial 
Cl ,  1  1; ••• ;5f 

6i  =  inherent  error  in  measurement  of  depth  or  pitch, 
i  =  1,  . .  • ,  5 ; 

ai  =  roundoff  error  after  addition  i; 
mi  =  roundoff  error  after  multiplication  i; 
a  =  roundoff  error  after  subtraction; 
fi  =  roundoff  error  in  gain  F(I,J); 

C6  =  normalizing  factor  to  convert  the  pitch  value  from 
a  16-bit  integer  word  to  a  compatible  value  with 
respect  to  depth. 

With 

|mantissa|  of  Y(l),  Y(2),  F(I,J),  Cl  to  C6  <  1.0 

and 

|ai|,  |  Ai  |  ,  |mi  |  ,  |fi|/  |a|  <_  2  16 
the  two  bounds  for  the  absolute  errors  become 
max | De ( 1) |  <  20  •  2-16  +  Ye(l)  +  Ye(2) 

<  1.25  •  2~12  +  Ye (1)  +  Ye(2) 

and  max|Ye(2)  |  _<  29  •  2  . 
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Assuming  for  Ye(l)  the  same  expression  and  combining  these 
inequalities  yields 

max|De(l)  |  £  1.25  •  2-12  +  2  •  29  •  2-16 
_<  0.001190  =  0.068  degrees. 

The  maximum  computational  error  due  to  the  16  bit  floating 
point  arithmetic  of  0.068  degrees  is  below  the  accuracy  of 
the  measurement  of  the  plane  deflection  (assumed  as  ±  0.5 
degrees).  The  error  in  D(I)  is  dominated  by  the  errors  in 
the  measurements  of  plane  deflections  and  therefore  the 
16  bit  floating  point  arithmetic  is  sufficiently  accurate 
to  carry  out  the  control  functions. 

3.  Utility  Routines 

For  debugging  and  display  purposes ,  some  additional 
routines  exist. 

a.  Routine  INTREAD  reads  four  hexadecimal  digits 
from  the  CRT  in  ASCII-code  and  converts  it  to  binary  code. 

b.  Routine  INTRNT  accomplishes  the  opposite  of 
routine  INTREAD. 

c.  Routine  FXPN TREAD  reads  a  fixed  point  decimal 
number  in  the  range  of  the  floating  point  format  from  the 
CRT  in  ASCII  code  and  converts  it  to  floating  point  format. 

d.  Routine  FXPNTPRINT  does  the  opposite  of  routine 
FXPNTREAD . 

e.  Routine  FLTPRNT  displays  a  number  in  floating 
point  format. 
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IV.  CONCLUSION 


This  investigation  has  shown  the  general  feasibility 
of  using  a  microcomputer  as  a  controller  for  depth  and  pitch 
of  a  submarine.  The  computational  speed  and  accuracy  of  the 
microcomputer  are  sufficient  to  permit  the  control  of  a 
submarine. 

In  addition  it  seems  possible  to  incorporate  the  control 
of  the  rudder  and  the  trim  in  the  controller.  Also,  some 
different  approaches  should  be  investigated,  i.e.,  the  use 
of  a  second  microcomputer  to  solve  the  Riccati  equation  in 
real  time  as  well  as  to  estimate  the  rates  of  change  by  a 
Kalman  filter. 

A  further  interesting  project  would  be  the  implementation 
of  the  submarine  simulation  on  a  microcomputer.  The  micro¬ 
computer  based  simulation  could  be  developed  for  use  in  a 
training  facility  to  exercise  manual  control  of  depth  and 
pitch  of  a  submarine. 
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APPENDIX  A 


EQUATIONS  OF  MOTION,  LINEARIZATION  OF  EQUATIONS 
OF  MOTION  AND  CROSS  REFERENCE  OF  NOTATION 

Al.  EQUATIONS  OF  MOTION 

The  following  set  of  equations  are  referred  to  a  body 
fixed  system  of  axes  which  are  coincident  with  the  prin¬ 
cipal  axes  of  inertia  of  the  body.  The  origin  of  this 
axis-system  is  located  at  the  assumed  center  of  mass  of 
the  body 

ation  of  Motion  Along  the  Body  Axis  System  x-Axis 


i  -  vr  +  wq)  ~  \  ^ 


1  q 2  +  X  1  r2  +  X 
n  rr  - 


*  u2  +  X  1  v2  +  X 
w 


ww 


vvn' 


'  (n'  -  1)  v2 


sin  0 


81 


nation  of  Motion  Along  the  Body  Axis  System  y-Axis 


v  -  wp 


+  ur,=  [y;  'r+Y.  -p] 


+  2l*  [Ypq  'Pq+  Yp|p|  ’Plpl] 

+  f^  [Yv  '*+Ywp  ,WP+  Yv|r|f|l(v2  +  w2)ll  hi] 
+  I*3  [Yr  ur  +  Y  |  r  |  5  r  'u|r|Sr+Yp  'up] 


+  2*3  Y  rn1  '  <n’  ’  1>ur 


+  |i2  [Y*  'u2+Yv  'uv+Yv|v|  '  v|(v2  +  w2)^|]- 


+  f*VY6r  '  5r 


+  l^Y6rn-  '  -D5r 


+  l^Yvn’  '  <n’ 


+  fi2Yv|v|n'  '  (n'  -  1)  v|(v2  +  w2)2  | 


+  ■?  1 2  Y  1  wv# 
£■  wv 


+  -f  X2  (F  )  v2  +  w2 

2  y  vs  — o — 


+  EW  |  sin  (j)  cos  0 


(- w)  Sin  cot 


#  Multiplied  by 


u 

U 


for  large  angles  of 
attack  near  -90° 
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Ojrq  Q-|(nj  °|(M 


Lon  of  Motion  Along  the  Body  Axis  System  z-Axis 


uq  +  vp)  =  -|  i4  Z  ^  1  q 


+  fA*  [Zrr 

1  r2  #  +  Z  ^  1  rp#J 

Note  1 

+  -&X3  [Zw 

1  w  +  Z  1  vr#  +  Z  1  vp  +  AZ  1  vp#l 

vr  vp  r  vp  ^  J 

+  Ii3  [Z  q  ' 

'  Uq+  Z  |q|6s  ,Ul^6S  +  ZW|q|’,i,l(v2  +  W2)ii  |q'] 

+  j  *3  Z  qn*  *  ^  ^  uq 

+  £*2  [Z  *  ’u2+Z  w  'uw  +  Zw|w|  ’WI(V2  +w2^l] 

+  £*’  f.Z|w|'«M  +  Zww'  lw  (v2+  w2)^|  +  Zvv'  v2#  ] 

+  !^2  [z  6s  '  5s  +  zab  ’6b] 

+  li2  [Z  wn*  ’  (n'  -  i)uw  +  Z  w|wln,  '  (n*  -  1)  w|  (v2  +  w2  )*  |  ] 

+  f*VZ6s»'  '  (n'  -  1)6s 

+  J?  X2  (F  )  V2  +  w2  V  Sin  cot 
2  z'vs  - - 

+  cos  0  cos  (j)  #  Multiplied  by 

u 

U 


for  large  angles  of 
attack  near  -90° 

Note  1 

when  not  multiplied  by 
u  add  to  Z  * 

TJ  VP 
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ation  of  Motion  About  the  Body  Axis  System  x-Axis 


+  (IZ-Iy)qr  =  pS 

\k  . 

1  p 

'  p  +  K 

qr 

'  qr  +  K  .  ' 

*  +  Kp|p| 

'  p|p|] 

+ 

[k 

L  p 

'  up  +  K  r 

'  ur  +  K  .  1 

V 

v  +  K  ' 

wp 

wp] 

*-§‘‘ 

[K*  ■ 

u2  +  K  ' 

V 

uv  +  Kv|v| 

1  v|(v2  +  w2 

♦l<* 

K 

vw 

1  vw 

+  -Si3  u2K  .  1  S 

2  6r  r 


+  Bz  g  sin  (j)  cos  0 
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i>n  of  Motion  About  the  Body  Axis  System  y-Axis 


Note  1 

j.  -  !z)  rp  =  •§  Is  [M.  '  q  +  M  rr  '  r2#  +  M  rp  '  rp  +  AM  rp  '  rp#] 

+  I l*  [M  q  '  u<*  +  M  |q|6s  '  u|q|6s  +  M  |w|q  '  |(v2+  w2^|q] 

+  -fz4  fM  .  'w+M  '  vr#  +  M  '  vp#l 
2  L  w  vr  vp  J 

+  f^qn'  *  (n*  ’  1)UC1 


+  i3  j^M  #  1  u2  +  M  ^  1  uw  +  M  w|  *w|(v2+  w3)2  |  J 

+  -?i3  Fm,  ,'u|w|+M  1  |w  ( v2  +  w2)2  I  +  M  '  v2  #  ”1 

2  L  |w|  11  ww  1  ’  1  vv  j 

+  f*3  u2  [M6s  'Ss  +  M6b  ’Sb] 

+  I"43  M  wn'  '  (n‘  -  !)  UW 

+  |-43Mw|w|n'  '  <n'  -  1)  w|(v2  +  w2)^  | 


+  f*3u2M6sn-  ’(n--  l)6s 


+  Bz  g  sin  0 


SW  .  xfc.  cos  0  cos  (j) 


§  Multiply  by  jj  for 
large  angles  ot  attack 


near  -90 


Note  1 

when  not  multiplied 
by  u  add  to  M  ' 

7  U  rP 
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uation  of  Motion  About  the  Body  Axis  System  z-Axis 


-Ix)pq  = 

[Nf 

’  r  +  N  '  pq  +  N  .  'pi 

pq  p  F  J 

+  | £ 1 

[Nr 

'  ur  +  N  |r|6r  ’  “|r|6r  +  N  |v|r  ’  |(v2  +  w2)^|r] 

+  il*  1 

N 

P 

*  up  +  N  .  '  v  +  N  *  wpl 

v  wp  r  J 

+  -2  X4  N  , 
2  rn1 

1  (n*  -  1)  ur 

*  i‘“  1 

>* 

•  u2+  Nv  'UV+NV|V|  '  v|(v2  +  w2)^  I J 

+  -|  X3  u2  N 

6r  '  6r 

+  fi3  u2N6rn-  '(n'-l)6r 

+  |i3  N  vn,  •  (n*  -  1)  uv 

+  |i3Nv|v|n,  •  (n1  -  1)  v|(v2  +  w2)^| 

+  -?  X3  N  •  wv# 

L  wv 

+  SW.  xti  cos  9  sin<f>  #  Multiply  by  -  for  large 

angles  of  atta<Hc  near  -90° 
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AUXILIARY  EQUATIONS 


p  +  i>  sin0 

(q  -  $cos0  sin<f)  )  /  cos({) 


=  (  r  +  0  sin<j)  )  /  cos  0  cos<J) 

=  ucos0cosij)  +  v  (  sinfy  sin  0  cos  ip  -  cos<j)  sini/) ) 
+  w  (  sin(|)  sin  ip  +  cos(j)  sin  0  cos  ip ) 

=  u  cos  0  sinj)  +  v  (cos  4>  cos  ip  +  sin(j)  sin  0  sin  ip  ) 
+  w  (  cos  ())  sin  0  sin  $  -  sin<t>  cos  ip  ) 

=  -  u  sin  0  +  v  cos  0  sincj)  +  w  cos  0  cos  4> 

=  (  u2  +  v2  +  w2  )2 


a i'  +  aa'  n'  +  a3'  n'2  J 

when  kx  <  n1 

§  **u2| 

[  V  +  V  a'  +  V  n'2  ] 

when  kg  <  n*  < 

f  ^u2 

[  cx'  +  Ca'  n'  +  c3'  n'2  ] 

j  when  kg  <  n*  <  kg 

[  dx'  +  ^'n'  +  d3'  n'2  _ 

when  n*  <  kg 

»  t  * 

aL»  a2>  a3 
i-  t  » 
bl,  b2>  b3 
*  *  » 
ci,  c2>  c3 

i  t  * 

dl ,  d2,  d^ 


Sets  of  non-dimensional  coefficients  used  in  the  pro¬ 
pulsion  equation  above.  The  set  which  will  be  in  effect 
at  any  time  during  a  simulated  maneuver  will  depend  on  the 
value  of  n'  and  the  numbers  k-^k 2*^3 ' 
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NOMENCLATURE 


All  symbols  used  in  the  equations  of  motion  and  in  the 
auxiliary  equations  and  relationships  which  appear  in  this 
report  are  defined  below.  Any  dimensions  involved  will  be 
consistent  with  the  foot-pound-second  system  of  units.  All 
angles  are  in  degrees. 


SYMBOL 


B 


SL 

u 


u 


V 


w 


DEFINITION 

A  dot  over  any  symbol  signi¬ 
fies  differentiation  with  res¬ 
pect  to  time. 

Buoyancy  force  which  is  posi¬ 
tive  upwards. 

Mass  of  the  submarine  including 
the  water  in  the  free  flooding 
spaces . 

Overall  length  of  the  submarine 

Linear  velocity  of  origin  of 
body  axes  relative  to  an  earth- 
fixed  axis  system. 

Component  of  U  along  the  body 
x-axis . 

Component  of  U  along  the  body 
y-axis . 

Component  of  U  along  the  body 
z-axis . 
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uc  Command  speed:  A  steady  value 

of  u  for  a  given  propeller  rpm 
whenc^ , ^  and  control  surface 
angles  are  zero.  Sign  changes 
with  propeller  reversal. 

x  Longitudinal  axis  of  the  body 

fixed  coordinate  axis  system. 

y  Transverse  axis  of  the  body 

fixed  coordinate  axis  system. 

z  Vertical  axis  of  the  body  fixed 

coordinate  axis  system. 

xQ  Distance  along  the  xQ  axis  of  an 

earth-fixed  axis  system. 

yQ  Distance  along  the  yQ  axis  of  an 

earth-fixed  axis  system. 

zQ  Distance  along  the  zq  axis  of  an 

earth-fixed  axis  system. 


p 

Component 

of  angular 

velocity 

about  the 

body  fixed 

x-axis . 

q 

Comp  onen  t 

of  angular 

velocity 

about  the 

body  fixed 

y-axis . 

r 

Componen  t 

of  angular 

velocity 

about  the 

body  fixed 

z-axis . 

ZB 

The  z  coo 

rdinate  of 

the  center 

of  buoyan 

ce  (CB)  of 

the  subma- 

rine  . 
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Angle  of  attack. 


(3 

Sb 


Sr 


Angle  of  drift. 

Deflection  of  bowplane  (  or 
sailplane  ) 

Deflection  of  rudder. 


SS  Deflection  of  sternplane, 

»  . 

n  The  ratio  uc/u. 


0 

$ 

<t> 

e 


CP 

t 


Pitch  angle. 

Yaw  angle. 

Roll  angle. 

Mass  density  of  sea  water. 

Weight  of  water  blown  from  a 
particular  ballast  tank  ident¬ 
ified  by  the  integer  assigned 
to  the  index  i. 

Angular  velocity. 

Time . 


x 


ti 


Location  along  the  body  x-axis 
of  the  center  of  mass  of  th£  i 

ballast  tank  when  this  tank  is 
filled  with  sea  water. 


th 
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<Fx> 


p 


I 


X 


Propulsion  force  (see  auxi¬ 
liary  equations  and  relation¬ 
ships)  . 

Moment  of  inertia  of  a  sub¬ 
marine  about  the  x-axis. 


I 


y 


Moment  of  inertia  of  a  sub¬ 
marine  about  the  y-axis. 


I 


z 


Moment  of  inertia  of  a  sub¬ 
marine  about  the  z-axis. 


Kp'.  V>  Kp|p|'>Kqr' 

Kr'  Kr'  Kv'  K*', 

Kv>  Kv|v|/»  ^vw  *  KSr/ 


Non-dimensional  constants  each 
of  which  is  assigned  to  a  parti¬ 
cular  force  term  in  the  equation 
of  motion  about  the  body  x-axis. 


I  ’  ,  M  '  ,  &M  '  ,  M  '  M,  |  ' 

rr  rp  rp  q  |  q  |  <5 s 

,  M.  1  ,  M  f  ,  M  1  ,  M  *  ,  M^  *  , 
r  w  *  vr  vp  qn  * 

I w I  *  M I W I  *  Mww  *  MVV  *  M6  s 

Mwn 1  *  Mw|w|n*  ’  M6sn’ 


Non-dimensional  constants  each 
of  which  is  assigned  to  a  parti¬ 
cular  force  terra  in  the  equation 
of  motion  about  the  body  y-axis. 


( 
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) ' , 

N 

1  ,  N.  1 

'  ,  N  1 

pq 

P 

r 

i'  , 

N.  * 

,  N  1 

1  ,  N 

V 

wp 

rn 

i|v 

» 

1  • 

N.  » 
6  rn 

W 

»  M  *  M  ' 


vn  v 


I  v  I  n ' 


Non-dimensional  constants  each 
of  which  is  assigned  to  a  parti¬ 
cular  force  term  in  the  equation 
of  motion  about  the  body  z-axis. 


X  1 

qq 


rP 


x. 1 

u 


wq 


Y  *  X  X  1  X  fX  1  Non-dimensional  constants  each 

uu  vv  ww  *  6  r5  r  *  6  s<5  s  * 

of  which  is  assigned  to  a  parti- 

X  1  X  *  X  ,  1  X  *  cular  force  term  in  the  equation 

o  bo  b  vvn  wwn  6s6snf  * 

of  motion  along  the  body  x-axis. 

Y  t 

A6  r6  rn  1 


Y.  :  ,  Y.  1 
r  P 


pq 


pIpI 


wp 


v | r |  ’  1 r  ’  1 1 r | 5  r ' ’  Yp'>  Yrn’’> 


'  Y  1  Y 
’  r  '  1 

Y*'  •  Y  '  >  Y  I  I  '  .  '  ,  Y.  ,  '  , 

*  v  v|v|’  5 r  6  rn 


Yvn  '  '  ’  Yv | v | n ' ' ’  Ywv’’  <Fy>vs 


Non-dimensional  constants  each 
of  which  is  assigned  to  a  parti¬ 
cular  force  term  in  the  equation 
of  motion  along  the  body  y-axis. 


2. ' ,  Z  1 
q  rr 


Zrp’>  V>  zvr> 


vp 


^Zvp  ’  Zq  ’  Z  j  q  I  6  s  ’  Zw|q| 

z-.. v. z z  i„r>  zi„r* 


z  ' ,  z  ' ,  z  '  ,  z,,  ' ,  z  ,  '  , 

ww  vv  os  6  b  wn 


Zw|w|n”  Z6sn '  ’  (Vts 


Non-dimensional  constants  each 
of  which  is  assigned  to  a  parti¬ 
cular  force  term  in  the  equation 
of  motion  along  the  body  z— axis. 
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A2.  EQUATIONS  OF  MOTION  SOLVED  FOR  u,  v,  w,  p,  q,  r 


The  equations  were  separated,  the  terms  involving  u, 
v,  w,  p,  q  and  r  on  the  left  side;  the  remaining  terms  on 
the  right  side.  The  equations  containing  Ix,  1^  and  Iz 
were  divided  by  V*  and  the  equations  containing  m  were 
divided  by  9? ,  where  Z  is  the  length  of  the  submarine, 
p,  the  density  of  the  fluid,  was  set  to  2.0.  This  yields 
a  system  of  equations  of  the  form 


(1)  A  •  DOT  =  Z 


where 


JSL  -  X!  0  0  0  0  0 

A3  n 


0  -Sy  -  Y!  -SY'  0  0  — AYI 

^  v  P  r 


DOT 


u 

V 

w 


p 

q 


r 


93 


and  Z  consists  of  the  remaining  linear  and  nonlinear  terms 

on  the  right  side.  Premultiplying  both  sides  of  (1)  by 

A  ^  solves  the  system  of  equations  for  u,  v,  w,  p,  q,  r. 

The  terms  and  -i  will  be  replaced  by  m*  and  I*  respectively. 
I6  JT 

In  contrast  to  the  notation  in  part  A1  of  Appendix  A,  the 
notation  for  the  speed  u  will  be  replaced  by  ua  to  distinguish 
it  from  the  notation  for  the  controls  u. 

A3.  LINEARIZATION  OF  EQUATION  OF  MOTION 

To  use  the  equations  of  motion  in  the  controller  the 
nonlinear  equations  were  linearized  under  the  following 
assumptions : 

i)  u  =  u  yields  u*  =  1,  all  terms  involving  (n'  -1) 
c  a  « 

vanish. 

ii)  The  linearized  model  is  in  trim,  therefore, 

\  wixti  cos  0  cos  =  0 

\  VT  cos  0  cos  =  0  . 

iii)  ua  is  a  measured  input, 

iv)  Depth  and  pitch  are  the  only  parameters  to  be 
controlled. 

The  remaining  nonlinear  terms  are  dropped.  Left  with  the 
equations  involving  w  and  q  the  problem  is  reduced  to  two 
degrees  of  freedom,  the  resulting  equations  are  shown  on 
the  following  page. 
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The  constant  additive  last  term  in  the  set  of  linearized 
equations  will  be  omitted  in  future  use  of  these  equations 
because  of  its  insignificant  contribution.  Let 


DET  =  (m*  -Z')  (I'  -M!)  -  Z'M' 
w  y  q  q  w 


and 


fm'  -  Z! 

w 

-HZ! 

q 

-I.-M.  HZ' 

INV  = 

1 

1 

< 

I'  -  M! 

y  q. 

_  1 
DET 

V*  m' -  z4 

then  (2)  becomes 


w 

’  vy* 

(m-  +Z^)ua 

w 

=  INV 

q 

_  uaMw/A 

uaMq/H  _ 

q 

+  INV 


UaZ5b/t 


u2aM6s/i.2  u2M&h/l2 


_  b_ 
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A4.  CROSS  REFERENCE  OF  NOTATION  IN  COMPUTER  PROGRAMS 

AND  IN  MATHEMATICAL  MODEL 

The  following  correspondence  between  the  notations  in 
the  mathematical  model  and  the  computer  program  has  been 
used: 

1.  Subscripted  variables  like  I  or  Z  i  i  are  written 

x  w  jq  | 

IX  and  ZW1Q1  where  the  1 1 s  denote  absolute  value 
signs . 

2.  Notations  which  are  obvious  from  their  use  in  the 
programs  are  not  explained  in  the  following  list, 
e.g.  IXY  =  IX  -  IY,  or  LC2  =  LC**2. 

3.  Some  of  the  notations  in  subroutines  can  be  traced 
back  through  the  argument  list. 

The  following' list  is  a  correspondence  between  some 
notations  and  explains  the  notation  where  necessary. 


Symbol  in 


Computer  Mathematical 
Program_ Model 


Description ,  if  necessary 


A(I,J) 


A 


Matrix  of  linear  coefficients  for  u,  v 


w,  p,  q,  r  of  equations  of  motion 


I,J  =  1  to  6 

AINV(I,J)  A-1  Inverse  of  A,  I,J  =  1  to  6 


3(1,  J) 


B 


Matrix  B  from  state  equatin,  I,J  =  1  to  4 


Cl 


Cl 


Element  of  weighting  matrix  S,  weight 


for  stern  plane 
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C2 


C1I 


C2I 


D(l) 

D  (2) 

CONO 
CONI 
CON  2 
CON  3 

DB 

DBER 


DBINC 


DBORD 


C2  Element  of  weighting  matrix  S,  weight 

for  fairwater  plane 

1/Cl  Inverse  of  Cl  combined  with  normalizing 
factor  for  D ( 1)  to  get  16-bit  integer 
representation  in  PLM 

1/C2  Inverse  of  C2  combined  with  normalizing 
factor  for  D(2)  to  get  16-bit  integer 
representation  in  PLM-program 

6g  Control  input  for  stern  plane  in  radians 

(used  in  controller) 

6^  Control  input  for  fairwater  plane  in 

radians  (used  in  controller) 

Normalizing  factor  for  0  in  PLM-program 
Limit  for  error  in  pitch  in  PLM-program 
Limit  for  error  in  depth  in  PLM-program 
Normalizing  factor  for  0  (1/CONO) 
in  PLM-program 

6^  Defelction  of  fairwater  plane  (used  in 

simulation  of  submarine)  in  radians 
Difference  between  actual  deflection  of 
fairwater  plane  and  control  input  for 
fairwater  plane 

Increment  or  decrement  of  deflection 
of  fairwater  plane  for  next  call  to 
routine  PLGEN 

Same  as  D(2) ,  but  in  16-bit  integer 
format  of  PLM-program 
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DBDOT 

DBGRA 

DEPTH 


DET 
DOT (1) 
DOT (2) 
DOT (3) 
DOT  (4) 
DOT (5) 
DOT (6) 
DPTH(I) 
DR 

DRINC 

DRER 

DRORD 

DS 

DSER 

DSDOT 

DSGRA 

DSINC 

DSORD 

F(I,J) 


Rate  of  change  in  movement  of  fairwater 
plane 

Same  as  DB,  but  in  degrees 
ZQ  Depth  of  submarine  (units  in  DSL/360 

program:  feet  -  units  in  PLM-program: 

inches) 

det  A  Determinant  of  A 
u 

v 

w 

P 

q 

r 


Last  five  values  of  DEPTH,  I  =  1  to  5 

<$r  Deflection  of  rudder  in  radians 

Same  as  DBINC,  but  for  rudder 
Same  as  DBER,  but  for  rudder 
Ordered  deflection  of  rudder 
Deflection  of  stern  plane  (used  in 
simulation  of  submarine)  in  radians 
Same  as  DBER,  but  for  stern  plane 
Rate  of  movement  of  the  stern  plane 
Same  as  DS ,  but  in  degrees 
Same  as  DBINC,  but  for  stern  plane 
Same  as  D(l) ,  but  in  16-bit  integer 
format  of  PLM-program 

F  Product  of  matrices  -K  and  B ,  I  =  1  to  4 , 


J  =  1,2 


99 


FF(IU,IW1, 

IW2 / I , J) 

Matrix  containing  all  values  of  F  for 

different  speeds  and  different  weights 

w^  and  w2 

G(I,  J) 

G 

Matrix  G  from  state  equation,  I,J  =  1  to 

IU 

Speed  UA  rounded  to  an  integer 

IW1 

Subscript  to  select  different  weights 

in  array  FF 

IW2 

Subscript  to  select  different  weights 

W2  in  array  FF 

K(I,J) 

K 

Gain  matrix,  I  =  1  to  4 

KIJ 

Same  as  K(I,J),  only  different  notation 

KDIG 

K 

Derivative  of  K(I,J)  as  a  function  of 

time,  I,J  =  1  to  4 

PERR 

Error  in  pitch  in  16-bit  integer  format 

of  PLM-program 

PIDOT 

0 

PITCH 

0 

PIGRA 

PITCH  in  degrees 

PORD 

r3 

Ordered  pitch 

PTCH(I) 

Last  five  values  of  PITCH,  I  =  1  to  5 

Q  (I  /  J) 

Q 

Weighting  matrix  for  error  in  pitch 

and  error  in  depth 

RODOT 

i 

ROLL 

S 

s 

Weighting  matrix  for  control  inputs 

SI 

s”1 

Inverse  of  S 
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UA 


UC 

W1 

W2 

Y(l) 

Y(2) 

Y  (3) 

Y  (4) 
YADOT 
YAW 
Z(I) 


ZERR 


ZORD 


u,u 


a 


W1 

w2 


zo  rl 


w 

Q"r2 


q 


rl 


Speed  along  body  x-axis  (notation  u  is 
used  in  simulation  of  submarine,  notation 
ua  is  used  in  controller  to  distinguish 
between  the  control  input  u  and  speed 
Command  speed 

Weight  for  error  in  depth,  element  of 
matrix  Q 

Weight  for  error  in  pitch,  element  of 
matrix  Q 

Output  variable  of  state  equation, 
error  in  depth 

Output  variable  of  state  equation 
Output  variable  of  state  equation, 
error  in  pitch 

Output  variable  of  state  equation 

Vector  consisting  of  all  nonlinear  terms 
involving  u,  v,  w,  p,  q,  r  and  remaining 
terms 

Error  in  depth  in  16-bit  integer 
format  of  PLM-program 
Ordered  depth 
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APPENDIX  B 


COMPUTATIONAL  PROCEDURE  FOR  THE  ERROR  PROPAGATION 
(shown  for  D(l)  and  Y(2)) 


The  relative  errors  for  addition  and  multiplication  are 
defined  as  [Ref.  8]. 

1.  Addition:  =  JL  +  JL  f* 

x+y  x+y  x  x+y  y 

e  e  e 

2.  Multiplication: 

x-y  x  y 

A  process  graph  is  a  graphical  representation  in  which  the 
arithmetic  operations  are  carried  out.  In  the  case  of 
addition  and  multiplication  the  scheme  of  labeling  is  the 
following: 
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To  obtain  the  total  error  leaving  a  node,  the  error  values 
entering  a  node  are  multiplied  with  the  branch  labels  and 
added  together.  To  this  sum  is  then  added  the  roundoff 
error  of  the  particular  operation. 

The  relative  error  Ye (2)  in  computing  the  rate  of 
changes  in  di  (see  also  process  graphs  at  the  end  of 
Appendix  B) . 


I  Ye  (2)  , 

I  Y(2)  I  1 


Al+61) 


Cl-dl 

Cl-dl+C2-d2 


Cl -dl+C2d2 
Cl-dl+C2d2+C4d4 


Cldl+C2d2+C4d4  .  Cl • dl+C2 • d2+C4 • d4+C5d5 
’ Cldl+C2d2+C4d4+C5d5  *  Cldl+C2d2+C4d4-C3d3 

Y  (I) 


C2  •  d2  ^  ,  A  „  ,  C4*d4 

+  (A2+62)  ■  -y  ^  j ^  +  (A4+54)  — y-fjJ 


,Af_lP,_v  C5-d5  ,  C3-d3 

+  CA5+65)  y  (ij"  +  (A3+S3)  Y(l) 


,  Cl-dl  , 
+  rnl  y'(x')’  + 


C2-d2 
Y  (I ) 


+  m4 


C4-d4 
Y  (I) 


+  m5 


C5-d5 

Y  (I ) 


+  m3 


C3-d3 

Y(I) 


Cldl+C2d2 

Y(I) 


a2 


Cl • dl+C2 • d2+C4 • d4 

Y(I) 


+  a3 


Cl • dl+C2 • d2+C4 • d4+C5 • d5 

YU) 


a 
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where 


Y ( I)  =  rate  of  change  in  depth  (Y(2))  and  in  pitch  ( Y ( 4) ) ; 

Ai  =  roundoff  error  in  Coefficient  for  Legendre  polynomial 


ci,  i  =  1,2, ... ,5; 


<5i  =  inherent  error  in  measurements  i,  i  =  1,2,.  ..,5; 

ai  =  roundoff  error  after  addition  i? 

mi  =  roundoff  error  after  multiplication  i? 


=  roundoff  error  after  subtraction; 


a 


di  =  sampled  data  (pitch, depth) ,  i  =  1,2,. ..,5. 

Multiplying  both  sides  with  Y(2)  and  rearranging  the  terms 
yields  the  following  expression  for  the  absolute  error 

|  Ye  (2)  |  £  |  Cldl  ( Al+<51+ml+al+a2+a3+a) 

+  C2d2  ( A2+<$2+m2+al+a2+a3+cr) 

+  C4d4  ( A4+<54+m4+a2+a3+a) 

+  C5d5  (A5+cr5+m5+a3+a)  +  C3d3 ( A3+63+m3-a) 

The  relative  error  De(I)  in  computing  the  ordered  plane 
deflection  is 


C6  *  Y  ( 3)  F  (4 , 1) 

D  (I) 


+  f 5+Ye (4) ) 


F  ( 3 , 1)  Y  ( 4 ) 

D  (I) 


+  (f 3+Ye (2) ) 


F  ( 2 , 1)  Y  ( 2) 
D(I) 


+  (f  l+<$2) 


F  ( 1 , 1)  Y  ( 1) 

dTTJ 


+  (ml+m2) 


C6-Y(3)F(4,1)  ,  F  ( 3 , 1)  Y  (4) 
- dTT) - +  m3  - Dll) - + 
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F  (2 , 1)  Y  (2)  ,  _c  F  (1 , 1)  Y  (1) 

D  (I )  +  ^  dTTI 


+  ctl 


C6-Y(3)F(4,1)+Y(4)F(3,1) 

D(I) 


,  C6*Y(3)F(4,1)  +Y  (4)F(3,1)+Y  (2)F(2,1)  , 

+  a J  D(x)  -  +  a2 1 


where 

fi  =  roundoff  error  in  gains  F(I,J); 

D (I)  =  control  output  for  plane  deflection; 

C6  =  normalizing  factor  for  variables  representing 
pitch  or  error  in  pitch. 


Multiplying  both  sides  with  D(I)  and  rearranging  the  terms 
yields  the  following  expression  for  the  absolute  error 

| De  (I)  |  _<  |C6-Y(3)F(4,1)  ( A6+61+al+a2+a3+ml+m2) 

+  F(3,X)-Y(4) (f5+Ye(2)+m3+al+a2+a3) 

+  F (2 , 1) -Y(2) (f3+Ye (4) +m4+a2+a3) 

+  F ( 1 , 1) *Y (1) (fl+62+m5+a3) | 
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Y  (I)  =  (Clxdl)  +  (C2xd2 )  +  (C4xd4 )  +  (C5xd5)  -  (C3xd3) 


Y (2)  PROCESS  GRAPH 


Cldl+C2d2 

Cldl+C2d2+C4i 


D(l)  =  (C6*Y(3))  •F(4,l)+F(3,l) -Y(4)+F(2,l)  •Y(2)+F(l,l)  *Y(1) 


D(X)  PROCESS  GRAPH 
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